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INTRODUCTION 


Time accurate and relaxation procedures for the compress- 
ible inviscid fluid conservation- law-equations (Euler equations) 
are frequently less efficient than procedures developed for 
potential flow. This is so because with the Euler equations 
one has to deal with a system of nonlinear partial differential 
equations as opposed to a single nonlinear partial differential 
equation in potential flow. The consequence of a system of 
equations is that a spectrum of associated eigenvalues must be 
considered so the equations are usually stiffer. Of more 
import, with a system of equations, the spectrum of associated 
eigenvalues can appear (as in subsonic flow) with both positive 
and negative real parts. This last restriction severely limits 
the choice of allowable difference operators and iterative or 
time dependent procedures. 

The purpose of this study was to develop and test relax- 
ation algorithms for the Euler equations which are made 
possible by use of flux vector similarity splittings. In this 
approach the flux vectors are split such that their associated 
Jacobian matrices have either all positive or all negative real 
parts. The split vectors can then be specially differenced 
using dissipative forward or backward operators. A major part 
of the study was devoted to improving the iterative convergence 
rate of such relaxation algorithms, using the multi-grid 
approach . 


BACKGROUND 


A Flux Vector Split Scheme 

The Euler equations can be written in cons ervative-vec tor 
form as 


8.q + 3 S + 9 F + 9 „G — 0 


X 


( 2 . 1 ) 


where 


-V 
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f -1 
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pu^+p 

pv 

> u = 

puv 
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u(e+p) 

• » 


etc . 


In at least one direction we can expect a convection velocity 
to be less than the local sound speed. Given this condition 
only a restricted mjmber of stable difference operators can 
be used. The problem is easily seen in one dimension if we 
assume subsonic flow, u < c where c is the sound speed. Then 


3tq + 3^E = 0 


( 2 . 2 ) 


By local linearization, since E = E(q) 

® + [If] (5 - 5o) 

=2 +A(q-q), A is the Jacobian matrix, 

o o ^ o 

Peculiar to the Euler equations is the relation E = Aq 
(the equations are first degree homogeneous) so 

2 A^q (2.3) 
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Thus, equation (2.2) can be locally approximated as 


^t^ + ^x^o^ ^ ^ (2.4) 

so that equation (2,4) becomes the linearized governing 
equation where is locally constant. 

By similarity transform 


A = S A S"^ or S"^A S = A 


where 


A 


"u 0 0 ■ 

0 u+c 0 

_0 0 u-c_ 


(2.5) 


are the eigenvalues of A, and S has as its columns the 
eigenvectors of A. 

Let U = q, then multiply equation (2.4) by and 

find using equation (2.5) and freezing 

+3 (A U) 
t o ^ 


or since A^ is diagonal 


9^U. + 3 (A.) U. 
t 1 x^ I'^O 1 


( 2 . 6 ) 


Because for subsonic flow A 2 = u+c > 0 and A 3 = u-c < 0, the 
operator 9 cannot be replaced by only a single difference 
operator unless that operator is a pure central difference. 
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All other operators will lead to exponential growth. However, 
first order central difference operators have problems. They 
tend to be nondissipative and they cannot be used with simple 
Euler explicit time differencing (i.e., point Jacobi relax- 
ation) . For these reasons it is difficult to devise good 
relaxation algorithms for the Euler equations using central 
differencing . 

Given equation (2.6), one would have no difficulty 
devising stable dissipative difference operators. If > 0 
let d be approximated by a backward differencing, if A. < 0, 
let 9 be approximated with a forward operator. Unfortunately 
it is impractical to try to diagonalize the Euler equations 
for this application. For the nonconservative equations, 
people have successfully differenced each term based on 
intuition, but this is an inappropriate way to proceed. 

A rational way to split the flux terms for special 
differencing is based on similarity transforms and the 
validity of local linearization. As noted previously 

f = A E [If] 

= SAS"'q 

there is no approximation here. 

The plus-minus (or positive-negative) splitting is as 
follows : 


A = + A" 

where 2A"^ = A + |A|, that is (if u > 0) 
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A"*" = 


= 


u 0 

0 u+c 

0 0 

u 0 

0 u+c 

0 0 


and a" = A - A"^. 


0 

0 

0 

0 

0 

U“C 


u < c 


u > c 


Then 


= SA'^' S ' q + SA~ S"^ q 

— + — — + 

= A q + A q 

++ ~ 

= + E 


(2.7) 


The flux vectors F and ^ can be split in the same way. For 
better clarification E is written out (2-D) 


= (p/y) 


(y-Dxt + .5(xt + xt) 

(y-l)uxt + .5u(xt + + -5c(xt - X^) 

(y-l)vX^ + .5v(xt + X'J’) 

. 5 (y-1) (u^+v^) X^ + . 25 (u^+v^+2cu+2c^/ (y-1) ) X^ 
+ .25(u^+v^ - 2cu + 2c^/(y-l))X^ 


(2.8) 


^ _ U+ 1 U I ^ + _ (u+c ) + I u+c 

Ai- 2 »A3- 2 > 


_ (u-c)+|u-c 

^4 0 


The vector E is obtained by replacing X"*~ by X”, while 
F"^ and F use X7=(v±|vl)/2, etc. 
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The advantage of the new splitting is that the eigenvalues 
associated with the Jacobian of E are alvjays positive while 
those of E are always negative. Consequently for the Euler 
equations split as 

3^_q + 9^(E^ + fi') +9 f~) +9 (G^ + G“) = 0 (2.9) 

t X y z 

one can easily pick stable dissipative spatial difference 
operators. For example, a stable point Jacobi relaxation 
scheme is given by 

-^n+1 ->n / c B , r B B 7 ^ , 

q = q - co (6 E+o F+o G + 

^ ^ ^ X y z 

+ 6 ^ 2" + 6 ^ F" + 6 ^ G“)^ (2.10) 

X y z 

where 

6 ^ u. = (3u. - 4u. + u. ^)/( 2 Ax^) 

X j ^ J j -1 3-2' 

6^ u^ = (-3Uj + ^^j+1 “ u^^_2)/(2Ax^) 

A stable, easily inverted implicit scheme is given by 

[I + h(5® A + 6 ® B"*" + 6 ® C'*')] [I + h(5^ A‘ + 6 ^ B" 

L y z -* X y 

+ C-)] - q*') = - h(&l e+ + r 

+ S® + 6 ^ E' + 6 ^ F” + S')” (2.12) 

z X y z ^ 
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The implicit algorithm was a chief motivation for 
+ and - splitting. The left hand side is factored into pure 
block lower and pure block upper triangular matrices and 
these are trivial to invert by simple sweeps (like SOR) . 
Compared to inverting block tridiagonal matrices in x, y, 
and z, the + and - inversion cost is very small. The differ- 
encing also does not require artificial viscosity. The dis- 
advantages of the + and - implicit scheme are that twice as 
many Jacobian matrices have to be formed (although A, B, and C 
can be formed all together making up some of the loss) , and 
no satisfactory way has been devised to time accurately treat 
viscous terms in upper and lower implicit factorizations. 

A more thorough discussion of the "upwind” flux vector 
split scheme is given in reference along with a review of 

previous related work. Additional new work in this area is 
given in references [ 2 ] through [sj . 

As an aside, it is noted that the fact that the Euler 
equations are homogeneous of degree one has been used in 
formulating the splitting scheme. For systems of conservation 
laws in which the flux vectors are not homogeneous, one can 
still proceed as before if one solves for a delta variable. 
Consider, for example, the system 

3^u + 9^f + 9yg = 0 (2.13) 

where f = f (u) , g = g(u) , but u f and u g . 

Let u^ be a nearby solution of equation (2.13), then 
expanding u, f, and g about u^. 
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® t [^0 + 

+ ®y[So + ®o('^-“oil = ° + O(u-u^)^ 

where A = ■— , ^ ^ fu * ^o satisfies (2.13), 

equation (2.14) is rewritten as 


(2.14) 


'tU + ^ ° 


(2.15a) 


and 


U = u-u 

o 


(2.15b) 


Equation (2.15) can now be split as were the Euler equations 


+ 


+ 


xijith A - A + A and B = B + B assuming that A and B 
can both be diagonalized. Note that because U = u-u^ is 
the new dependent variable, the scheme (2.15) will only be 
first order accurate in time if A is set to A and U is 

■n4-l n ^ n4-1 

set to u^ - u“. If A i^ set to A^ and U is set to 
(3u^ - 4u^ + u^ ■^)/2, second order accuracy can result. 

For steady state problems one can use u-u^ = u^-Uj^^, etc. 
if, for example, u^ is known at j = 1 and u^ satisfies the 
steady form of (2.13) . This technique has not been tested 
in steady cases. 


It is stressed that no computational experience has been 
gained with split forms of equation (2.15). While + and - 
split forms of equation (2.15) are not as computationally 
efficient as (2.9), they apparently will not suffer from the 

9 E*! I 

annoyance that ^ A^ etc. as discussed in Ref. 1. 
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Multigrid 


A major emphasis of this study was to apply the multigrid 
relaxation methodology ^7j - ||lo] to the flux vector split form 
of the differenced Euler equations. The multigrid procedure 
has been notably successful in treating elliptic second order 
partial differential equations, and recently has been applied 
to the transonic potential equations [9]> [11] , . Consider- 

ably less experience is available with first order systems of 
partial differential equations. 

The multigrid method has been described in various 
publications, and it is clear from these papers that the 
method is still in a state of change. The work of Brandtl , 
who is perhaps multigrid's most distinguished proponent, 
treats multigrid as a smoothing process in which coarse grids 
are used to damp low frequency components of a solution and 
fine grids are used to damp high frequency content. A recent 
paper by Jameson [jll] used multigrid with an alternating 
direction algorithm in which the various grids effectively 
supply a sequence of relaxation parameters to damp the 
spectrum of frequencies associated with a given solution. 

In another recent work McCormick [13] has analysed multigrid 
procedures as a means of supplying efficient conditioning 
matrices to a relaxation procedure. The concept of high and 
low frequency content is again observed in this analysis. 

With the exception of some very simple model problems, 
the multigrid procedure is too complex to treat by analysis. 

(The analysis should account for the interpolation process, 
boundary conditions, and so on.) Consequently, one must 
proceed using simplified analysis and trial and error. 

The danger here is that it may be difficult to identify 
the crucial elements of a successful multigrid procedure. 
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The multigrid procedure used in the current vjork is 
essentially that which has been referred to as the full 
approximation. The general forms used are described 

below. However, there are various strategies for sequencing 
through and interpolating from grid level to grid level. 

These have not been faithfully followed in the results to 
be presented later. 

Let G^, ...., G^, G^ be grids with ever 

increasing finess such that G 2 is half the spacing of G^, 
is half the spacing of G 2 , etc. Using multigrid sweeps 
the variables are updated on the mth grid as 

C" = QS + Hm (2.16) 

where the residuals are used directly rather than their local 
linearizations . 


Here q represents the best fine grid (i.e., G^) estimate 
of the exact difference equation solution 


RrXq) = 0 = 6^ E"*" + 6^ E" + 6'^ F"*" + 6^ F 
X y y y 

The variable is the update variable. At the end of each 
raultigrid sweep q is replaced by Qm' while in changing from 
grid to grid the discrepancy is updated: 


^+1 


m-f-l 

•M 


+ I 


m+l 

"m 


(Qm - 
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The matrix H is a conditioning matrix which is formed as 

part of the chosen relaxation scheme, while is the 

9 R nH"l ^ 

Jacobian matrix, — . The operator means to interpolate 

m-1 


The operator I 


m 


9q ■ m 

data from coarse grid m to finer grid m+l . 

means to average (or inject) data from fine grid m to coarser 

grid m-1. 

As an alternate to equation (2.16), the update formula 




QS + - An.(q)q + iS 


Rw^q)] 


(2.17) 


has also been used. This last form, which is more costly 
to form, has the nonlinearity frozen over a multigrid sweep. 

The multigrid procedures defined above essentially 
duplicate the full approximation scheme of South and Brandt. 
Motivation for the various steps in this procedure can be 
found in reference [9] or in Appendix A. In Appendix A, the 
basic steps of the multigrid process are derived using simple 
expansion and interpolation ideas. 
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ONE- DIMENSIONAL FLOW STUDY 


Nttmerical Schemes 


If one is careful to avoid direct- inversion fully 
implicit schemes, it is possible to test relaxation algorithms 
for the Euler equations in one dimension. The one dimensional 
nozzle flow equations are given by 


3 F + H = 0 

X 

where 


■p ■ 


-pu 


-0 - 

pu 

, F = A 

pu^+p 

, H = -p 

^x 

_e 


_u(e+p) 


_0 _ 


(3.1) 


(3.2) 


3 A 

Here A is the area ratio and A is r — . In one dimension, 

X o X 

pressure is related to the other variables via 


p = (y- 1) (e - %pu^) 


(3.3) 


Using flux vector splitting, (3.1) is differenced as 


s’’ F'*' + 6^ F' + H = 0 = R(q) 

X X 


where 


V- - PA 

^ " 2 ? 


2(Y“1) Xj + Ag + X 3 

2(y-l)uA“ + AJ(u+c) + A~(u-c) 


L (Y-1) AJ u^ + 


Aj(u+c)^ Ag(u-c) 


-H 


+ W 


(3.4) 


(3.5) 


where 


W = 


(3-Y)(aT + aT)c 
2Ty^T) 


xf = !i-4- 


u 


etc . 


(3.6) 
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I 


and 


X 


F = 


3F. - 4F.., + F.., 
2Ax 


F = 


-3F. + 4F.^, - F.^, 

2Ax 


(3.7) 


Two iterative schemes have been tested for equation (3.4). 
The first is point Jacobi 


^n+l ^n /2Ax^ ,3F 

Ij >^3q 


,+ 


3F“v“^j^n 


3q 


(3.8) 


where — and — are the Jacobian matrices and q is defined 
3q 3q ^ 

in (3.2). The second is an approximate factorization (AF) 
scheme 


+ 

/T , xb 3F ^ , xf F ^ /-n+1 -n^ 

(I + 00 6 ) ( I + 00 0 ) (q ~ q ) 

^ x3q^^ xq ^^ ^ ' 


= - 00 R 


n 


(3.9) 


Note that the second scheme could be solved without approximate 
factorization in the block pentadiagonal form 


(I + oo 6 


b 3F 

X- 9q 


,+ 


, j~f 3F~s //-n+1 

+ ('5 - q ) = 


- u Rj (3.10) 


with 00 >> 1. Such a process cannot be efficiently repeated 
in multidimensional, however, so only the factored form 
(3.9) is considered. The form (3.9) is representative of 
multi -dimens ions since essentially the same type of lower 
and upper factored forms are obtained in two or three 
dimensions . 


13 


For comparison purposes a centrally differenced form 
of (3.1) was also considered. In this case the flux vector 
F was split into a convection term and a pressure or sound 
speed term, i.e. 

F = F + F (3.11) 

u c 

/N /N /-V /N 

with F^ = uq and ^ splitting was 

arbitrarily introduced for the purpose explained below. 
Introducing (3.11) into (3.1) and central differencing 
gives 

6 F + 6 F + H = 0 = R (q) (3.12) 

XU X c c^^^ 


where 


6 F = 

X 


F - .1 - F. T 

= j+J- j-i 


2 Ax 


This set of difference equations was then solved via the 
implicit approximately factored scheme 

fl + (0 6 a1^ + (e + + lVAq'^|J7A] [l + m 6 a” 

L. XU A.C. 

+ (e + h^ + |VAq^|)VA] (q^"*"^ - q^ ) 

= -uR(q") + (eVA + + [VAq'^DVAq^ (3.13) 

^ J 

where 


u 


’ 0 

1 

0 


u 

0 

0 ' 

-u^ 

2u 

0 

or A^ = 

0 

u 

0 

ue/p 

e/p 

u 


0 

0 

u 


(3.14a) 
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0 0 
Aq = (Y-1) u^/2 -u 

(u^-ue/p) (e/p-|- u^) 

and ntomerical dissipation terms are added to the equations 
with h~Ax, e = 0(.06), and Vq^ = q^ - etc. This 

way of treating the dissipation terms was developed over 
several years, but see also Appendix B. 

Note that equation (3.12) could have been ^^itten as 
one factor and then easily inverted via a block tridiagonal 
solver. With w this amounts to a direct inversion of 

the linearized equations. The splitting equation (3.11) 
allows use of the approximately factored form, equation (3.13) 
which is the type of factorization required in multi-dimensions. 
Thus, it is suspected that information learned from study 
of the scheme equation (3.13) will apply to the multi- 
dimensional case. One other point --the matrices A and A 

3Fc ^ ^ 

are not the Jacobians and — , because use of the 

dq dq ' 

Jacobians leads to numerical instability. However, the 

4 J via similarity 
eigenvalues of 
c. The eigen- 
values of A^ are u, u, u, while those of A^ are 0, c, -c. 

One Dimensional Results 

One dimensional nozzle flow was selected as a test. 

The nozzle was chosen with area ratio 

A = 1 - .8 X (1-x) 0 £ X 1 


matrices A and A , developed in reference 1 

dF ^ ^ *- 

to , are stable. It is remarked tha t the 

dFu ^ x: SFq / y — 1 

are u, u, u, those of are 0, 
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x^hich has a throat to exit ratio of 0 . 8 . The entrance area 
ratio at X = 0 is equal to the exit area ratio at x = 1. 

The incoming Mach number must reach a value of 0.55332 for 
the throat to reach sonic conditions. In all one dimensional 
test cases, boundary data at both entrance and exit is fully 
specified. This overspecification of boundary data is 
expected to aid the iterative rate of convergence. 

Subsonic flox«7 convergence rates for plus and minus flux 
vector split schemes, equations (3.8) and (3.9), are indicated 
in figure (3.1). Here the maximum residual is plotted as a 
function of fine grid iteration number. The £2 norm of the 
residuals was also monitored, and this showed a similar, 
although somewhat more monotone, trend. A sixty five point 
grid was used. On this grid the subsonic flow results are 
indistinguishable from the exact solution on a typical plot. 

It was hoped that the multigrid scheme would make even 
the point Jacobi scheme converge quickly. As indicated in 
figure (3.1) this was not the case. Of course, some modifi- 
cation of the current multigrid scheme may lead to a quite 
different result. In this and all of the one dimensional 
calculations the multigrid scheme (2.16) was used. 

The point Jacobi scheme was quickly dropped as it 
requires as much computational work per step as the AF 
scheme, equation (3.9). The AF scheme converges well by 
itself. With use of multigrid, however, the rate of convergence 
is improved by a factor of 4, (see figure (3.1)). Of course, 
the multigrid scheme requires more work per point, but this 
result is promising. It is noted that in the multigrid 
approach that sweeps from coarsest grid to finest grid 
proved more efficient than sweeping down from the finest 
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grid and then back up from the coarsest to the finest grid. 
Linear interpolation was used to go from a coarse to a fine 
grid. Injection (not averaging) was used to go from a fine 
to a coarse grid. By injection is meant that the fine grid 
result is used on the coarse grid, as the coarse grid points 
occupy the same position as the fine grid points. 

Numerical results for supercritical shocked flow are 
indicated in figure (3.2) versus the exact solution. The 
small "glitch" at the sonic line is due to a discontinuous 
change of the A and A eigenvalues. This effect would be 
more pronounced, but smoothing has been added by simply 
replacing A“ by 

_+ + 

A “A ± e 


where e is small. This effectively adds 


e (6 


b 

X 


6^)q 


E (Ax) " 


6 


xxxx 


q 


to the equations, where 6 is the five point centered 

difference approximation to the fourth derivative. 


The rate of convergence for the AF scheme with and without 
multigrid acceleration is indicated by figure (3.3). Again a 
sixty five point grid is used. While the multigrid scheme ends 
with a much faster rate of convergence, it does not reach a 
converged state any faster than the standard algorithm. 


The central differencing scheme was also run for this flow. 
Numerical results are shown in figure (3.4), while the 
convergence rate is indicated via figure (3.5). Here the 
multigrid procedure gives a superior rate of convergence to 
the conventional scheme, and overall, it nearly matches the 
rate achieved with plus and minus flux vector splitting. 
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The use of the dissipation terms are crucial to the success 
of this method. While this way of treating the dissipation 
terms was derived over some years, their treatment is 
motivated from the plus and minus flux vector split teirms 
as indicated in Appendix B. 
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TWO DIMENSIONAL FLOW STUDY 


Test Problem and Numerical Scheme 


As a test problem in two dimensions, the flow about a 
nonlifting biconvex airfoil was used. Thin airfoil boundary 
conditions were specified so as to simplify the imposition 
of the tangency condition. As a result it was unnecessary 
to introduce transformations to fit the body surface. 
However, simple clustering transformations of the form 


5 = ?(x) 

n = n(y) 


(4.1) 


were used so as to stretch the far field boundaries far 
from the airfoil. The transforms also provide good grid 
resolution near the body, especially at the leading and 
trailing edges of the biconvex profile. 

With introduction of the C,n variables the governing 
equations take on the form 


+ E‘) + 3^1 + F') = 0 (4.2) 

where E“ was given by equation (2.8) and F“ are similar 
(see reference |[l] to obtain the actual elements) . Note 
that E and F remain the original Cartesian form of the flux 
vectors, only the independent variables are transformed. 

A multigrid relaxation algorithm is developed for 
equation (4.2) using the approximate factorization scheme 
(2.9) and the full approximation multigrid method (2.17). 

The precise algorithm used is 





= “ 

m m 


qS> 

(4.3) 


where 






and A , B"*”, E"*”, etc. are evaluated using q. 

The parameter h^ is a relaxation parameter which varies 
from point to point as 


= Jh , 
m m ’ 


h = 0(0-2) 

m m 


(4.4) 


A more conventional scaling would use 


\ = \ “ 0 <^> <^- 5 ) 

and perhaps (4.5) or alternation of (4.4) with (4.5) would be 
preferred to (4.4). 


All of the far field boundary conditions for the 
nonlifting biconvex airfoil are set to free stream values. 
A stretched grid is used and these boundaries are placed 
3 chord lengths in front of and behind the airfoil. The 
upper boundary y = is set at 8 chords away from the 

y = 0 plane. Along the y = 0 line the vertical velocity v 
is obtained via 


V 


U, 


dy 

dx 


(4.6) 
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where ^ = 2t(1-x), 0 <. x £ 1, and ^ “ 0 elsewhere. 
All other variables along y = 0 are obtained from |-3L = q. 

°y 

In the multigrid scheme this implies that 




m 


^^j,l 


(Qm - 


m 

M 




(4.7) 


Results 

Subsonic flow results about a six percent thick parabolic 
arc airfoil are shown in figure 4.1. The Mach number is 
0.817. These results, which use the thin airfoil boundary 
condition, are comparable to those obtained in potential 
schemes. The numerical solution shows slight as 3 mimetry 
about the midchord and this may be taken as a measure of 
error in the solution. A stretched grid with 65 points 
in X and 33 points in y is used. 

Convergence rates for this case without using multigrid 
acceleration are indicated in figure 4.2 for various values 
of h^. The best multigrid results obtained to date are 
indicated in figure 4.3. The multigrid procedure is perhaps 
four times faster and gives a fully converged solution in 
less than 50 fine grid iterations. It is remarked that 
linear interpolation and injection (not averaging) is used 
in the multigrid procedure. It is also noted that the grid 
spacing does vary rapidly. In particular, ^^max*^ ^^min ~ 36.06 
and Ay_ /Ay.„* = 44.2 so grid aspect ratios are of 0(40). 

Supercritical flow results about a 127o thick profile 
are shown in figure 4.4 for various smoothing ratios, i.e. 

± X ± I A I 

\ = 1 — L ± e , e a constant 

2 


Here the Mach number is 0.8. 



Convergence rates obtained without using multigrid are 
given in figure 4.5, the best multigrid results appear in 
figure 4.6. The multigrid result appears to be about twice 
as fast as the nonmultigrid result and it requires about 
120 fine grid iterations using the grid of 65 x 33. Perhaps 
the multigrid procedure is less effective for transonic 
flow because of the extensive discontinuous switching that 
occurs in the split flux vector formulation. 



CONTINUOUS DERIVATIVE PLUS -MINUS SPLIT FLUX VECTORS 


As the results in both one and two dimensions illustrate, 
small ’’glitches” appear at the sonic line when the plus -minus 
flux vector splitting is used. This is because the flux 
vectors depend on the |a|, as indicated via equation (2.8) 
or (3.6). As A changes sign, |A| has a discontinuous first 
derivative. Consequently, the flux vectors also have a 
discontinuous first derivative when any eigenvalue changes 
sign. This occurs at the sonic line and causes numerical 
inaccuracy there, a glitch. 

In reference [U it was noted that the A^ should be 
redefined as 


_ A_j|; — \_X_\_ + 3 ^ 3 positive (5.1) 
2 

where ”3 is a small positive number which smoothly approaches 
zero as |A| increases.” In this way A“ could be kept smooth, 
however, no simple formula for 3 was proposed. 


The appearance of discontinuous first derivatives in 
the definition of the flux vectors prompted van Leer of 
ICASE to postulate new flux vectors. In one dimension he 
proposed (unpublished) 


(u - c) < 0 





pc (M+1) ^ /4 

f"*" = 

fa+ 

- 

pc^M+l)^ [(y-1)M+2]/[(4y)] 


^3 _ 


Y^ [(ft)Vf|]/[2(Y^-l)] 

F' = I 

- F 

+ 



(5.2a) 


(5.2b) 
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and 


u - c > 0 F"*" = F , F' = 0 (5.3) 

Here 



pcM 


pu 


F = 

pc^ (M^ + y”^) 

= 

PU^ 

+ p 


PC^[(MV2) + (y-1)“^]m 


u(e 

+p) 


These vectors have functional variation much like those 
proposed by Steger and Warming ^ij . However, they have a 
continuous first derivatives at u = c and at u = 0. Analysis 
by van Leer verifies their stability as do our own numerical 
tests. On the nozzle problem the van Leer vectors give the 
excellent steady state result shown in figure (5.2) which 
has a smooth transition across the sonic line. The second 
order upwind special difference operators, (3.7), were used. 

The success of van Lber ' s scheme prompted another 
search for a clean way to define 3 in equation (5.1). 

This was ultimate ly achieved by simply replacing the 
|X| with ^ e As shown in figure (5.1), this amoimts 
to fitting a hyperbola to the |x| versus A curve. 

i 

Thus A are redefined as 


A” 


A ± Va^ + 

2 


(5.5) 
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r 


and 3 of equation (5.1) is indeed positive and equal to 

(5.6) 

+ 

Use of (5.5) to define A" also leads to split plus minus 
flux vectors with continuous derivatives. One dimensional 
nozzle flow results using this definition of X are indicated 
in figure (5.3). Here the variation at the sonic line is 
smooth, while compared to the previous result figure (3.2), 
the shock solution is essentially unchanged. 
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CONCLUDING REMARKS 


Relaxation algorithms were devised for the steady gas 
d 5 mamic equations in conservation- law- form using flux vector 
splitting and the multigrid technique. Overall, fairly good 
steady state convergence rates were obtained although the 
multigrid method was not made as effective as what appears 
to be possible. Considerable additional work remains to be 
undertaken. 

In closing it is noted that other promising approaches 
for relaxing the Euler equations are beginning to appear . 
These include the surrogate equation approach of Johnson [1.5 
As a final digression it is suggested that a time-accurate 
second degree wave equation formulation, as detailed in 
Appendix C, might offer another entertaining approach. 
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Fig. 4.1. Subcritical solution obtained with plus-minus 
flux splitting. 
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Fig. 4.4. Solution for flow about biconvex using 
plus-minus flux split scheme. 
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APPENDIX A 


Multigrid Methodology 

Denote the difference equations for the flxix split 
partial differential equation as 

R(q) = + 6 ^F" = 0 

X X y y 

where 


E = E(q), F = F(q) 

Define grids Gj, G 2 , G 3 , . . . , G^, . . . , where G^^ is the 

finest grid. As in matrix notation, the higher the value 
of m, the more elements or points. Let q be a solution to 
the difference equations, and let q be an approximate solution. 
Expand R(q) = 0 about q, then 

R(q) = R(q) + A(q-q) + O(q-q)^ (A.l) 

where A = is the Jacobian matrix of the entire set 

9q 

of difference equations. 

On the fine grid (A.l) is rewritten 


^M^‘^M”^M^ ^M^^^ ” ^ (A. 2) 

M-1 

where subscript M denotes the finest or Mth grid. Let 

be the averaging operator from M to the grid M-1, then 

M-1 

(A. 2) can be multiplied by to obtain 

^M ^M^^M ■ ^M^ "" ^M " ^M^ 
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r 


or approximately 


4 '^ ^M-l - V <A.3) 

M 

where interpolation operator from grid M-1 to M. 

Note that 


•M -pM-1 ^ ,M 

M-1 M M 


(A. 4) 


where is an approximate identity of rank corresponding 

with the Mth grid. 


To the same approximation one can let 


M-1 

M 


A. 


M 



^M-1 


9R 


M-1 


'^M-1 


(A. 5) 


Although this last step is unnecessary, it is useful for 
programming ease. With use of (A. 5), equation (A. 3) is 
rewritten 

^M-1 ^M ^*^M ■ ^M^ "" ^M ^^M “ (A. 6) 

which can be represented as 

Ae = f (A. 7) 

where A = coefficient matrix, e = discrepancy variable, 
and f = residual forcing function. 

Without making the substitution (A. 5), the only required 
approximation in the reduced system (A. 6) or (A. 7) is in 
equation (A. 4). Except for this approximation, a solution 
to the reduced system (A. 6) is a solution to (A. 2). 


43 



Rather than solve (A. 6) or (A, 7) directly, an iterative 
solution can be sought of the form 

- e’^ = H(Ae"-f) (A. 8) 

where H is the conditioning matrix defined by the iteration 
s cheme . 

Once an approximate iterative solution is obtained from 
(A. 6), call it the new estimate of is obtained 

as 

M 

" ^M-1 ^M-1 


where 

M-1 M-1 M M 

1 _1 

and Aj^ ^ is an iterative approximation to typically 

n-1 

n=0 

is the first few terms of the Neumman series . This overall 
process can be carried out over more than one grid, and with 
use of various iteration strategies, defines the multigrid 
technique . 

Note that the above discussion doesn't show that this 

process is efficient or even desirable. It simply defines 

V M M-1 

the mechanics. It seems clear, however, that 1) 1^^ 

should approximate I as closely as possible; and 
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2) that whatever matrix is used to approximate » it 

1 M 

should be as close as ^M-1 possible. Examples 

of the interpolation matrices follow. 

M+1 

Interpolation Operator 


Given data on the odd indexed points, as shown, 

000000000 

123456789 

use linear interpolation for even indexed points, that is 

Ui + U3 U3 + U5 

U 2 = , Uj^ = , etc. 

2 2 

In matrix form the equations are written for a nine point 
grid as 


u; 


' 2 


Ui 


Ui 

U2 


1 1 


U3 


(U2+U3 ) /2 

U3 


2 


Us 


U3 

U4 


1 1 


U7 


(U 3 +Us )/2 

Us 

= % 

2 


Ug 

= 

Us 

Ue 


1 1 


(US+U7) /2 

U7 


2 


U7 

Ub 


1 1 


(U7+Ug)/2 

,Ug_ 


2 


Ug 


In the multigrid notation the equations are represented as 

T-m+l 

m+1 m m 

1^ I T 

Use of other than linear interpolation defines a new matrix I 
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Averaging Operator 


Given data at all grid points , 

000000000 

123456789 


average the data for the interior odd-indexed points, 
that is 


u, = 


U2 + 2Ug + U5 


U 5 = 


+ 2U5 + Ug 

_ 


, etc . 


In matrix form the equations are written for a nine point grid 
4 


Ui 

U 3 

Us 

U 7 

Ug 


= % 


12 1 
12 1 
12 1 


4 


Ui 


r 

Ul 

U2 


(U 2 + 2 U 3 +U 4 ) / 4 

U 3 

= 

(U4+2Us+Ug ) /4 

U4 


(U6+2u7+Us) /4 

Us 


Ug 

U6 


L 

U7 


Us 


Ug 



In the multigrid notation the equations are represented 


as 


U = I , T U , T 
m m-fl m-f 1 
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If injection is used the matrix 
5 x 9 


has the form for the 




0 

0 10 
0 10 
0 10 
0 1 


Approximate Reduced Identit-' 


1 ^TTi ___ ^ 

m m-1 m-1 



8 0 0 0 0 

16 10 0 
i 0 1 6 1 0 

0 0 16 1 

0 0 0 0 8 
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APPENDIX B 


The relation of the upwind plus -minus split scheme to 
central differencing schemes is reviewed here following 
unpublished work of Ami Harten. We previously defined 

A+ = S s-' , A-' = s ^ s-' 

so 


_ S A S”^ L s A 
2 2 

= (A + |A|)/2 , = (F + |A|q)/2 (B . la) 

and 

A' = (A - lA|)/2 . F" = (F - |A|q)/2 (B.lb) 


Apply these into our equation (8.4) of TM 78605 
(reference Ed) 

[ • n A _L XI I ^ ^ 

I + (V a: , I + A A7 , I + V B T , I + A BT 1 i )lAq7^ 

l + 5^xj,k' Xj,k^ yj,k' yj,k' aI^j, 

= -(t- 4V ^ (<5^Ft , 1^ + 6^FT , I + 6^ct , I + 6^GT , ^ 

1 + " X j,k' X j,k' y j,k' y j,k' 


+(rjrj ) Aq";^ 


(B.2) 
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giving (pull % out in front) 


+ Vy|B| - Ay|B|)]Aq” = - | i^) («x<F + 

+ 6^(F -lAlq) + 5^(G-HB|q) + 6^(G -[Blq)}" + ^ Aq""'^ 

(B.3) 


Define the centered difference operators 


V + A E - E 

X X 


- 1 


2 Ax 


= 6 


V . A = - . g~^'± 21 - E-^ ^ -(Ax)-^AV) 
XX Ax \ \ / 


5^ + 6^ 

X X 


+E"^-4E"^+4E-E’^^ 


4 Ax 




X 


(B.4a) 

(B.4b) 

(B.4c) 


6^^ - 6^ = E'^-4E-4 6-4E + eIA ^2Ax) ^ (VA) ^ (B.4d) 

^ ^ 2Ax 

and substitute these into (B.3) 


[l+ (|||) («x^"+'5yB" - (2Ax)'‘(VA) |A| - (2Ay)'\vA) |B|)] Aq’’ 
= - (j^) (5 jj.f’^+7^g”+(4Ax) ■ ‘ (VA) ^ I A| q"+(4Ay) ' ' (VA) ^ | B | q"^) 

+ 14:|- 4q”'^ (B.5) 


Equation (B.5) represents a conventional central difference 
scheme with second and fourth order dissipation terms. 
Further simplification is possible. 
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Define the spectral radius p = l^max^ replace 

A I with p^, |B| with p^ . Then (B.5) is approximated as 

VAp, 


[ 


I + 0^ (6 + 6 b" - 

(1+^) ' X y 2 Ax 


)]Aq’" 


- (t4f) (« f" + 6 G’^ + 
^ 1 +^ ^ ^ X y 


2Ay 

, ('''') 'Pb .n 
4 Ax 4Ay 


q ) + ) ^q 


n-1 


(B.6) 


or in factored form 
9At / c- *n 


[ 


I + 


TTkT 


X 


VAp 

a 

2 Ax 


)][ 


0At 


VAp, 

n 


I+(T% --2Ay 


)] Aq" 


= RHS (B.6) 


(B.7) 


Equation (B.7) has the same type of numerical dissipation 
terms used in our conventional central difference scheme 
(c.f. Steger and Baily, AIAA J., March, 1980). 


The left-hand side of the scheme (B.2) can be factored 
into the product of two operators as 


[I 


+ h(V At , 1^ + V Bt , 1^)1 ri + h(A A. , 

^ X j ,k' y J L ^ X j ,k 

+ )] = RHS (B.2) 


where 


h = 


9At 

1 + e 


(B.8) 
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Substitution of (B.l) into the above 


[l + h(V^(^±^)+ Vy(?±i?i)] [l+h(A^(^^^) +Ay(^ 
= RHS (B.3) 


Replace |A| = p^I and |B| = p-j^I which degrades the ti 
accuracy to 0(At) . Then 


[l + hV 


(pi + A) hV, (p, I + B) 


+ hA 


X 2 

(B - P, I) 


+ 




][ 


I + 


hA^(A - pj) 


y 2 


] Aq" = RHS (B.6) 


^ )] Aq^ 
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APPENDIX C 

Using a "Lax-Wendrof f trick” one can write the Euler 
equations as a second order wave-equation. This leads to 
a time accurate formulation which for steady state problems 
is a particular case of Johnson's surrogate equation approach. 

The Euler equations are differentiated in time 

+3E+3F)=0 

giving 

(3. + 3 A + 3 B) 9+.q ~ 0 

t X y ^ t^ 

as 3^3^E = 3^3^.E = S^E^q^ = 3^Aq^ , etc. 

But 3 q = -3 E + 3 F so (C.2) is rewritten as 

(a similar substitute is used in the Lax Wendroff scheme) 

3^^q = (3 A + 3 B)(3 E + 3 F) (C.3) 

tt^ Vx y^^x y ' 

Expanding (C.3) gives the desired form 

3^^q = 3 A3 E + 3 A3 F + 3^B3 E + 3 B3 F (C.4) 

tt^ X X xy tx yy 

or making use of the fact that the equations are homogeneous 
of degree one 

(3^^ -3A3A-3A3B-3B3A-3B3 B) q = 0 (C.5) 
^tt XX xy yx yy 

The last form (C.5) shows the structure of the equation while 
the right hand side of (C.3) gives a particular case of 
Johnson's surrogate equation approach. 
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